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Abstract 

We present a first-principles study of 180° ferroelectric domain walls in tetrag- 
onal barium titanate. The theory is based on an effective Hamiltonian that 
has previously been determined from first-principles ultrasoft-pseudopotential 
calculations. Statistical properties are investigated using Monte Carlo sim- 
ulations. We compute the domain-wall energy, free energy, and thickness, 
analyze the behavior of the ferroelectric order parameter in the interior of the 
domain wall, and study its spatial fluctuations. An abrupt reversal of the 
polarization is found, unlike the gradual rotation typical of the ferromagnetic 
case. 
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The cubic perovskites are among the most important examples of ferroelectric materials^ 
Many undergo not just one, but a series, of structural phase transitions as the temperature 
is reduced. These transitions occur as a result of a delicate balance between long-range 
dipole-dipole interactions that favor the ferroelectric state, and short-range forces that favor 
the high-symmetry cubic perovskite phase. Because of the anomalously large Born effective 
charge of the atoms, the ferroelectric transitions in the perovskites are very sensitive to 
electrostatic boundary conditions.!'! As a consequence, domain structure plays an important 
role in the ferroelectric transitions, and a theoretical understanding of the domain walls is 
of great interest. 

Theoretical investigation of ferroelectric domain walls has been much less extensive than 
for their ferromagnetic counterparts. The strong coupling of ferroelectricity to structural and 
elastic properties is problematic. Previous theoretical investigations have concentrated on a 
phenomenological level of description, using Landau theory to study domain-wall thickness 
and energy.Bi Simple microscopic models such as local-field theory have also been used to 
identify the domain-wall structure and character.li Due to the limited experimental data 
available, this empirical work has tended to be qualitative and oversimplified, and has thus 
not been able to offer the accuracy needed for a deeper theoretical understanding. 

In this paper, we undertake a first-principles study of domain walls in BaTiOs. To our 
knowledge, this is the first such study in any ferroelectric material. Using an ab initio 
effective Hamiltonian developed previously to study the phase transitions of BaTiOs,! we 
set up Monte Carlo (MC) simulations to investigate the structure and energetics of 180° 
domain walls of (100) orientation. In particular, the energy, free energy, and thickness of 
the wall are calculated. We also analyze the behavior of the ferroelectric order parameter 
in the interior of the domain wall, and study the fluctuations in the domain-wall shape. 
Where we can compare with previous work, we find our results in general agreement with 
experiment aliHi3 and theoretical^^ reports. 

Because only low energy distortions are important to the structural properties, we work 
with an effective Hamiltonian written in terms of a reduced number of degrees of freedom.0 



The most important degrees of freedom included are the 3A^ "local-mode amplitudes" Uia 
for site i and Cartesian direction a. A "site" is a primitive unit cell centered on a Ti atom, 
and the "local mode" on this site consists of displacements of the given Ti atom, its 6 nearest 
oxygen neighbors, and its 8 nearest Ba neighbors, in such a way that a superposition of a 
uniform set of local-mode vectors Uj = e (independent of i) generates the soft zone-center 
ferroelectric mode polarized along e. We also include six degrees of freedom to represent 
homogeneous strain of the entire system, and 3A^ displacement local-mode amplitudes Via 
that serve to introduce inhomogeneous strains. We thus reduce the number of degrees of 
freedom per unit cell from fifteen to six, simplifying the expansion considerably. 

Since the ferroelectric transition involves only small structural distortions, we repre- 
sent the energy surface by a Taylor expansion around the high-symmetry cubic perovskite 
structure, including up to fourth-order anharmonic terms where appropriate. The energy 
consists of five parts: an on-site local-mode self-energy, a dipole-dipole interaction, a short- 
range interaction between local modes, an elastic energy, and a coupling between the elastic 
deformations and the local modes. The Hamiltonian is then specified by a set of expansion 
parameters, which are determined using highly accurate LDA calculations with Vanderbilt 
ultrasoft pseudopotentials.0 The details of the Hamiltonian, the first-principles calculations, 
and the values of the expansion parameters have been reported elsewhere.i This scheme has 
been successfully applied to single-domain BaTiOa to predict the phase transition sequence, 
transition temperatures, and other thermodynamic properties with good accuracy. 

The phase transition sequence for BaTiOa is cubic to tetragonal to orthorhombic to 
rhombohedral as temperature is reduced. We focus on the tetragonal phase, since it is the 
room-temperature phase, and the best studied experimentally. We adopt the convention 
that the polarization, and thus the tetragonal c-axis, are along z. In this phase, two kinds 
of domain walls, 90° and 180° walls, are possible.0 The notation refers to the angle between 
polarization vectors in adjacent domains. We choose the 180° domain wall for this study 
because of preliminary indications of a simple structure and narrow width.S Because it is 
energetically unfavorable to form domain walls carrying net bound charge, 180° domain 
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walls are restricted to lie parallel to the polarization. Earlier work has indicated that the 
180° domain wall of (100) orientation has much lower energy than for other, e.g. (110), 
orientations.il Thus, we focus on the (100) domain wall, which we take to lie in the y-z 
plane. 

Previous work has indicated that the width of the 180° domain wall is very narrow, of 
the same order of magnitude as the lattice constant .Sili For a very narrow domain wall, 
our choice of local mode (Ti-centered as opposed to Ba-centered) may introduce some bias. 
The point is that the sharpest domain wall that can be constructed is one for which the 
local-mode vectors Uj are constant except for a sudden sign reversal from one plane of sites 
to the next. For the Ti-centered choice of local modes, this represents a Ba-centered domain 
wall, for which the atomic displacements have odd symmetry across (and vanish on) the 
central Ba plane. Conversely, for a Ba-centered choice of local mode, the sharpest domain 
wall is Ti-centered, vanishing on a central plane of Ti atoms. In order to determine which 
of these scenarios is the more realistic, we constructed 4x1x1 supercells (containing 20 
atoms and two domain walls) corresponding to each of the above scenarios, using a mode 
amplitude taken from the average equilibrium structure of the MC simulations (very close to 
the experimental structure). We then performed LDA calculations to compare the energies 
of the two structures. We find the Ba-centered and Ti-centered walls constructed in this 
way have energies of 6.2 erg/cm^ and 62.0 erg/cm^, respectively. Thus, a sharp Ti-centered 
domain wall appears very unfavorable, and it is clearly best to use a Ti-centered local mode as 
we have done. We note that the effective Hamiltonian reproduces the energy of this sharpest 
Ti-centered wall to within 1% of the LDA result (not surprisingly, since configurations of a 
similar kind were included in the fitting§). 

We study the structure and energetics of the domain walls using Metropolis MC 
simulations.^ The degrees of freedom are the vectors and v. for each site z of a 4L x L x L 
supercell, and the six homogeneous strain components. As discussed below, the supercell 
is arranged to contain two domains, each roughly of size 2L x L x L, with domain walls 
normal to x. Periodic boundary conditions on the supercell are assumed. Since all energy 
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contributions (except for the dipole-dipole coupling) are local, we choose the single-flip MC 
algorithm. We make a trial move of variables at one site, check acceptance, make the change 
if accepted, and go on to the next site. One Monte Carlo sweep (MCS) constitutes one entire 
pass through the system in this manner. 

To generate a reasonable starting configuration, we equilibrate an L x L x L supercell at 
a high temperature (T >400 K) in the cubic phase, and then cool it down slowly, allowing it 
to relax for 20,000 MCS's at each temperature step. We stop the cooling when the tetragonal 
phase is reached, in which the polarization vector averaged over the simulation cell points 
along one Cartesian axis. (As reported in Ref. |^, this phase corresponds to the temperature 
range from 230-290 K in our calculation, while the actual experimental range is 278-403 
K.0) If the polarization is not along +5, we rotate the structure to make it so. We then 
copy the structure four times along the x-axis, with the polarization reversed to —z for two 
of them, as shown in Fig. |I|. 

The supercell is thus constructed to contain two periodic 180° domain walls perpendicular 
to the (100) direction. This structure is initially equilibrated for 2,000 MCS's, and then 
thermodynamic averages are constructed from runs of 40,000 MCS's. Since a supercell with 
domain walls has a higher energy than one with just a single domain, we have found that the 
domain walls sometimes spontaneously disappear during long simulations. To prevent this, 
we fix the 2;-components of the u vectors in the central two layers in each domain during the 
simulations. Since this structure should be very close to equilibrium and the constrained 
components are far from the domain wall, we think the effect on our results is negligible. 

Fig. 1^ shows a snapshot of the polarization vector components averaged over y-z layers, 
Ux, Uy, and function of x, for L = 10. Several qualitative features are immediately 

apparent. First, the sharp reversal of Uz indicates that the domain boundary is indeed very 
sharp, its width being on the order of a lattice constant. Second, the other components 
and Uy remain small throughout the whole supercell, and their random fluctuations do not 
appear to be correlated with the domain-wall position. (The qualitative difference between 
the fluctuations of Ux and Uy with x is an artifact of the averaging and of the presence of 
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strong longitudinal correlations.llj) Thus, we find that the domain boundary entails a simple 
reversal, rather than a rotation, of the ferroelectric order parameter. 

These behaviors are to be contrasted with the case of ferromagnetic domain walls, where 
the magnetization vector typically rotates gradually (on the atomic scale), keeping a roughly 
constant magnitude. This difference in behavior can probably be attributed largely to 
the much stronger strain coupling in the ferroelectric case. For our BaTiOa geometry, for 
example, the entire sample, including the interface, develops a tetragonal strain along z, 
imposed by the presence of domains polarized along ±z. This gives rise to a strong anisotropy 
which will tend to keep the ferroelectric order parameter from developing components along 
X or ?/ in the interface region. Thus, instead of rotating, the polarization simply decreases 
in magnitude and reverses as we pass through the domain wall. This absence of rotation 
of the polarization has been experimentally verified for the case of the 90° domain wall in 
BaTi03.lli 

We now turn to a quantitative analysis of our simulation results, focusing on the domain- 
wall width, smoothness, and energy. We first estimate the domain-wall thickness t as follows. 
For a string of sites along x at a given value of {y, z) and on a given MCS, we identify the 
pair of sites between which Uz changes sign. We then define t via the linear extrapolation 
t/a = 2u^^°^^ /Auz, where a is the lattice constant, u^p"^^ is the spontaneous polarization 
deep in a domain, and Auz is the change of Uz between the two interface sites. Finally, 
we average over {y, z) points and over MCS's to get an average value of t. The value of 
t estimated in this way is 1.4 unit cells, or 5.6 A. This is in reasonable agreement with 
empirical theoretical estimates of 6.7 A (Ref. §), and experiments which place an estimated 
upper bound of 50 A (Ref. H). 

To analyze the smoothness of the domain wall, we Fourier transform the polarization Uz 
as a function of the x coordinate for each {y, z) point, and retain only the first three terms 
in the expansion. This is an effective way to smooth the data while keeping the most useful 
information. The positions of the two domain walls in the supercell, denoted by Xi and X2, 
are identified with the values of x at which the Fourier-smoothed Uz changes sign. In this 
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way we obtain Xi{y, z, r) and X2{y, z, r), where r labels the MCS. 

In Fig. 1^, we show the probability distribution of Xi for L = 10 and for a run of 40,000 
MCS's at 260 K. The solid line is a histogram of the values of Xi{y, z, r), while the dashed 
line is a histogram of y-z planar averages Xi(r). A comparison of the two curves shows 
that the spatial fluctuations of the domain-wall position are much smaller than its ensemble 
fluctuations. From the solid line, we see that the Xi values have a typical standard deviation 
of between 1 and 2 lattice constants. Other runs indicate that this result is not very sensitive 
to system size. So, we can conclude that the domain walls are relatively smooth. We can 
further separate the contributions to these fluctuations coming from the y- and z-directions. 
It is found that the fluctuations along the 2;-direction (i.e., along the polar direction) are 
about 40% smaller than along the y-direction. The sign of this result was to be expected, 
since the shape of the domain wall should be such as to minimize the surface charge AP • n 
that develops on it. Here, AP is the change of the polarization vector across the domain 
wall, and n is the unit vector normal to the wall. 

Finally, we turn to an estimate of the domain-wall formation energy. Because of the 
periodic boundary conditions imposed on our system, there are no surfaces to give rise 
to a depolarization energy. Thus, the domain-wall energy can be calculated from the 
difference between the energy of the ALxLxL supercell with and without domain walls. This 
difference is small, but because the correlation time of the system (far from the transition) 
is quite short (20 MCS's), a sufficiently long simulation is capable of reducing the statistical 
errors in Ey^ to an acceptable level. The calculated domain-wall energies are shown in Table 
|. The reported values have a statistical uncertainty of about 4%. Simulations for two 
lattice sizes and temperatures are reported. We can see that our results are well converged 
with respect to system size. Because of the large increase of the correlation time near the 
transitions, it has proven difficult to give accurate values for Eyj at other temperatures. 

Our calculated value of Ew=16 erg/cm^ for the domain-wall energy is, however, probably 
not the proper quantity to compare with experimentally derived values. Instead, we should 
compute a free energy, E^, which includes entropic contributions from fluctuations of the 
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ferroelectric order parameter in the vicinity of the domain walL A glance at Fig. 0, which 
shows considerable fluctuations, suggests that such contributions are likely to be important. 

We have estimated the domain-wall free energies Fyj using an adiabatic switching tech- 
nique, as follows. First, we start with an equilibrated AL x L x L supercell containing two 
domain walls, and for which the ^-components of the u vectors in the central two layers in 
each domain are constrained to preset values, as before. We slowly reverse the values of the 
constraint variables in the center of one of the domains over the course of a 20,000-MCS 
simulation, making a small change in the constraint variables every 10 MCS's, and compute 
the total work done on the constraint variables. If the simulation succeeds in removing the 
two domain walls adiabatically, we can equate the work done to twice the domain-wall free 
energy Fy^. By comparing runs of from 20,000 to 30,000 MCS's, we find differences in com- 
puted values of only about 10%, which suggests that the switching is indeed adiabatic. 
The resulting computed values of F^, are about 4-5 erg/cm^, or about 3-4 times smaller than 
the values (and slightly smaller than the 6.2 erg/cm^ reported above for the energy of 
the ideal Ba-centered wall). 

Our result is consistent with previously published estimates of the "energy" of the (100) 
180° domain wall, although such previous values are rather scattered and inconclusive. Pre- 
vious experimental results of 10 erg/cm^ and 3 erg/cm^ (close to the phase transition) 
were given by Merzl^ and Fousek and Safrankova,0 respectively. On the theoretical side, 
Bulaevskiili reported a value of 10.5 erg/cm^ using a continuum Landau model, while 
Lawless! calculated an energy of 1.52 erg/cm^ based on a microscopic phenomenological 
model. (Since the above estimates involve use of empirical models fit to finite-temperature 
data, the "energy" values are probably best interpreted as free energies.) 

This investigation has opened several avenues for further studies. Using the same method 
described here, it should be possible to carry out a similar analysis for other types of domain 
walls, e.g., (110) 180° or 90° domain walls. Bigger MC simulations on larger systems would 
allow the computation of more precise values for the energy and thickness of the domain 
wall. 
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In summary, we have studied the properties of 180° domain walls in BaTiOs using a first- 
principles based approach, by applying Monte Carlo simulations to a microscopic effective 
Hamiltonian that was fitted to ab-initio total-energy calculations. The simulations were 
carried out in the middle of the temperature region of the tetragonal phase, relatively far 
from the C-T and T-0 transitions. We confirm that the domain walls are atomically thin, 
and that the order parameter does not rotate within the wall. We quantify the width, 
smoothness, and energetics of these domain walls. Our theoretical values of the wall width 
and free energy are in reasonable agreement with previously reported values, where available. 

This work was supported by the Office of Naval Research under contract number N00014- 
91-J-1184. 
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TABLES 

TABLE I. Calculated domain-wall energies Eyj and free energies Fyj as a function of simulation 

cell size L and temperature T. Statistical uncertainties are about 4% for Eyj and 10% for F^. 



T (K) L Ey, (erg/cm^) (erg/cm^) 

250 8 15.8 4.6 

260 8 17.1 4.0 

250 10 15.6 5.0 

260 10 17.0 4.4 
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FIGURES 



FIG. 1. Schematic illustration of the simulation supercell arrangement. 

FIG. 2. Snapshot of the y-z layer-averaged polarization- vector components Ux (dotted line), Uy 
(dashed line), and Uz (solid line), as a function of x/a (a is the lattice constant), for the 40 x 10 x 10 
lattice at 260 K. 

FIG. 3. Histograms of domain-wall positions for the 40 x 10 x 10 lattice at 260 K. Solid line, 
histogram of Xi{y, z,T)/a values (a is the lattice constant, r labels a MCS); dashed line, histogram 
of y-z planar-average values Xi(r)/a. 
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